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Projected Hartree-Fock theory (PHF) has a long history in quantum chemistry. PHF is liere 
understood as the variational determination of an A'"-electron broken symmetry Slater determinant 
that minimizes the energy of a projected state with the correct quantum numbers. The method 
was actively pursued for several decades but seems to have been abandoned. We here derive and 
implement a "variation after projection" PHF theory using techniques different from those previously 
employed in quantum chemistry. Our PHF methodology has modest mean-field computational 
cost, yields relatively simple expressions, can be applied to both collinear and non-coUinear spin 
cases, and can be used in conjunction with deliberate symmetry breaking and restoration of other 
molecular symmetries like complex conjugation and point group. We present several benchmark 
applications to dissociation curves and singlet-triplet energy splittings, showing that the resulting 
PHF wavefunctions are of high quality multireference character. We also provide numerical evidence 
that in the thermodynamic limit, the energy in PHF is not lower than that of broken-symmetry HF, 
a simple consequence of the lack of size consistency and extensivity of PHF. 
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I. INTRODUCTION 

In a recent paper^i we introduced a novel wavefunc- 
tion method known as projected quasiparticle theory 
(PQT). The fundamental idea of PQT is very simple. 
We begin with a Hartree-Fock-Bogoliubov determinant 
(thereby mixing states of difFererent particle number) 
and then restore particle number symmetry with the 
aid of projection operators in a self-consistent variation- 
after-projection (VAP) manner. Without breaking any 
other symmetries, this reduces to the standard antisym- 
metrized gcminal power (AGP) wavefunction. However, 
we can also deliberately break and restore spin symmetry, 
point group symmetry, and complex conjugation symme- 
try, all within the same general framework. In each case, 
this deliberate symmetry breaking and restoration allows 
us to construct fairly sophisticated multiconfigurational 
wavefunctions with mean-field computational scaling. 

Closely related to PQT, at least conceptually, is pro- 
jected Hartree-Fock (PHF), where particle number sym- 
metry is not broken and the projection operators act on 
an iV-electron determinant. The PHF acronym in quan- 
tum chemistry is usually associated with spin-projection 
in a projection-after-variation (PAV) approach^ where 
the deformed {i.e., symmetry broken) determinant \^) 
variationally minimizes the energy 
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and is then used to construct a projected state 
^1$) whose energy is given by 



($|pli?P|$) _ ($|i7P|$) 
($|Ptp|$) ~ (<I>|P|$) 



(2) 



the Hamiltonian. The PAV approach is computationally 
quite simple, but has several drawbacks. Near regions 
of spontaneous symmetry breaking, the PAV energy can 
be ill behaved.— Moreover, the projection is often carried 



out only approximately)^ 



without removing all spin 



Here, we have used the fact that the projection oper- 
ator P is Hermitian, idempotent, and commutes with 



contaminants. This may lead to a large deviation in the 
expectation value (S*^) for the projected state when the 
reference determinant has contaminants of many differ- 
ent spinsi^ Lastly, the wavefunction is not variationally 
optimized, complicating the evaluation of properties de- 
pending on its derivatives. 

These defects can be remedied by using a self- 
consistent VAP approach, wherein one obtains the de- 
formed determinant |$) by variationally minimizing the 
energy of Eqn. [2] with respect to |$), as first pro- 
posed by Lowdin in his extended Hartree-Fock method 
(EHF)i^ More often than not, EHF has been associ- 
ated with the use of a spin projection operator on a 
reference unrestricted determinant (the so-called spin- 
projected EHF—), and we will follow this usage. We 
should note that Goddard's GF method^ is equivalent 
to EHF, while in certain limits also the spin-coupled va- 
lence bond (SCVB) method^ corresponds to EHF. The 
use of projected wavefunctions in quantum chemistry is 
nowadays very limited, but they are ubiquitous in fields 
such as nuclear physics.—"— 

Several authors have proposed different ways to op- 
timize the EHF wavefunction in quantum chemistry. 
Mayer's derivation was based on the Brillouin theorem 
for the projected state.— His derivation relied heavily 
on the pairing theorem by LowdinJ^ Rosenberg and 
Martinoi^, and later Klimo and TinoiS., used a direct 
minimization of the energy functional. More recently, 
Byrman suggested an optimization technique based on 
the fact that the EHF wavefunction is recovered when 
a single spin wavefunction is used in the context of the 
spin-coupled valence bond methodi^ Our purpose in this 
paper is to show how to carry out a VAP optimization of 
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a PHF wavefunction in an alternative, computationally 
efficient, manner. In achieving this goal, we are guided 
by our previous work.— 

The derivation of the PQT equations relied on ex- 
pressing a general Hartree-Fock-Bogoliubov determinant 
1$) using its Thouless parametrization with respect to 
the bare vacuum |— ). In other words, we implicitly 

wrote 1$) = e'^j— ).— This parametrization requires 
($1— ) 7^ 0. In the case of PHF, however, the deformed 
determinant is a number eigenfunction with N electrons 
and is thus orthogonal to the bare vacuum. The PQT 
derivation cannot then be blindly applied. The practical 
consequence of this orthogonality is that one encounters 
indeterminancies in the PQT equations when using refer- 
ence determinants whose one-particle density matrix has 
occupation numbers equal to zero. 

While the PQT derivation cannot straightforwardly be 
extended to derive corresponding PHF equations, the 
underlying conceptual framework is similar. That is, 
we wish to express the energy of the projected state 
in terms of the one-particle density matrix of a broken 
symmetry reference determinant. Having done so, the 
energy can be minimized with respect to the reference 
determinant's density matrix, leading to an effective one- 
particle Hamiltonian to be variationally optimized. We 
thus follow the basic ideas used in PQT to derive the 
PHF equations from scratch. Our derivation and result- 
ing equations afford spin symmetry restoration of both 
coUinear and non-coUinear broken symmetry determi- 
nants in a straighforward manner, where coUinear (unre- 
stricted Hartree-Fock-type) determinants are eigenfunc- 
tions of Sz and non-collinear determinants (generalized 
Hartree-Fock-type) are not. Additionally, we can readily 
restore complex conjugation or point group symmetry as 
well, just as we did in PQT. All of this is accomplished 
in essence by choosing a more convenient representation 
of the projection than the standard Lowdin projection 
operators used historically. 

The use of non-collinear reference determinants in the 
optimization of projected Hartree-Fock states has been 
limited to the work by Lunell on the two-electron series.— 
The complex molecular orbital method of Hendekovioi^ 
is closely related to our complex conjugation projection 
but has seen little use. Unlike these approaches, how- 
ever, our implementation of PHF readily allows for the 
simultaneous restoration of many symmetries at once. 
Our method is also related to the VAMPIR method of 
Schmid et al.,— though the optimization of the reference 
determinant is carried out in a very different way. 

The properties of the projected Hartree-Fock wave- 
function have been studied extensively by many authors. 
We point the interested reader to the review by Mayer— 
on the spin-projected EHF method as a starting point 
to access the wide literature on the subject. In partic- 
ular, we mention that properties of the density matri- 
ces characterizing the standard EHF states (i.e., those 
based on coUinear determinants) have been derived by 
Harriman<2i Hardisson and Harrimaui^, Mestechkiui^ 



and Phillips and Schugi^ Simons and Harriman^^ de- 
rived properties of the point-group projected Hartree- 
Fock wavefunction. 

We should also point out some of the vices of PHF: 
the method is neither size consistent nor size extensive, 
as it has been noted before;^ Recall that a size consis- 
tent method is one where the dissociation limit of the 
system AB is equal to the energy of A plus the energy 
of i?; a size-extensive method is one where the energy 
is proportional to the number of particles. Our results 
in this paper confirm previous conclusion o^^'^'' that, in 
the thermodynamic limit, the PHF energy per particle 
reduces to that of the broken symmetry mean field. In 
other words, as the number of electrons becomes large 
enough, PHF has nothing to add over standard Hartree- 
Fock except that it retains good quantum numbers which 
the broken symmetry Hartree-Fock state loses. 

Section |TI] discusses the details of the projection we 
use before deriving the energy of our projected wave- 
function and the variational equations we solve to ob- 
tain it. In Sec. IIIII we present the results of applying 
our equations in the description of molecular dissocia- 
tion processes and in the computation of singlet-triplet 
splittings, before drawing a few conclusions in Sec. IIVI 
In the interest of brevity, we have deferred the detailed 
equations for our effective Hamiltonian to the appendix. 



II. THEORY 

A. Projection Operators 

The complicated nature of Mayer's EHF equations 
owes much to the form of the projection operator used. 
This form is due to Lowdin,— and writes 



P^=n 



l^s 



8^-1(1 + I) 
s{s + 1) ~ l{l + 1)' 



(3) 



Applying this operator to a wavefunction composed 
of multiple different spins just returns the component 
with spin s, and annihilates the components with other 
spins. While this spin operator does not project the z- 
component of spin, it can be generalized to do so. Given a 
simpler way of writing the projection operators and eval- 
uating the projected energy Eqr in terms of the deformed 
determinant |$), one would have a computationally sim- 
pler scheme whose results are identical to those of Mayer 
in the case of spin projection. 

Fortunately, such a simple representation of general 
projections existsJ^i^ In our previous paper we omit- 
ted some details about the form of the projection, a fact 
we intend to remedy here. We start our discussion by em- 
phasizing that by projected Hartree-Fock theory (PHF) 
we mean recovering a wavefunction with good quantum 
numbers from an intrinsically deformed HE state, and do 
not demand the use of actual projection operators in a 
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strict mathematical sense. In other words, while we re- 
quire that the wavefmiction be an eigenfunction of 
the relevant symmetry operators, we do not insist that 
P is either Hermitian or idempotent. 

If a single operator A can be associated with a constant 
of motion (that is, A commutes with the Hamiltonian) , 
one can write a projection operator as 

P^^y f d0e'*(^-^), (4) 

where A is the eigenvalue recovered and L is the corre- 
sponding volume of integration^^ We note that this form 
can be used to project onto eigenfunctions of Sz ■ 

If, on the other hand, there is a set of operators that 
commute with the Hamiltonian (but not necessarily with 
each other) and which form a group G = {g}, then one 
can recover the correct symmetries by diagonalizing the 
Hamiltonian in the basis formed by the elements of the 
group. In other words, 

|*)=^c,5l$), (5) 
a 

with the coefficients being determined by the solution to 
the corresponding eigenvalue problem. 

Alternatively, one can recover the coefficients by de- 
manding that the wavefunction j'l') has the desired sym- 
metries. In other words, one can work with projectors. 
In general, one can construct the so-called "transfer" 
opcrators22i2£ 

PL -^iY.^' (9)1. 9, (6) 

9 

where h is the order of the group G, Ij is the dimen- 
sion of the irreducible representation , and T^{g)xK. is 
the element in the A-th row and K-th column of the ma- 
trix associated with g in such irreducible representation. 
This operator yields zero unless the function on which it 
acts belongs to the K-th row of F^ . By using the great 
orthogonality theorem one can show that 

iPiJ = PL (8) 

It is easy to see that Pjl^ is indeed a projection operator in 
the mathematical sense: it is Hermitian and idempotent. 

Consider the action of P;^^ on a deformed wavefunction 
1$). It extracts the component of |<I>) which transforms 
as the K-th row of . It is, however, unphysical in the 
sense that 

PLm^PLi'^') (9) 

where |$') is a rotated wavefunction in the subspace of 
. In order to avoid this unphysical behavior, we use 
the linear combination 

I*) = E^-^aJ*)- (10) 



It is clear that this linear combination produces a wave- 
function I^P) which transforms as the A-th row of F-' , thus 
having A and j as good quantum numbers. 
Evaluating the energy of j'J') leads to 

^ ^ E..'<,c4mpLyHPij'f) 
_ E..><,c4mp^,,j'P) 

E.K'<'C.($|P;^;j$) ' 

where we have used the fact that P^^ commutes with the 
Hamiltonian. The coefficients are most conveniently 
determined by solving the generalized eigenvalue problem 

T.K'.c.=Ej2<'.c., (12) 

K K 

with hi,^ = ($|7?Pf,„|$) and = ($|Pf, 1$). Ob- 
serve that all projected wavefunctions of the same irre- 
ducible representation F-' are degenerate. 

We have used the above projections to restore spatial 
symmetry. On the other hand, we have preferred to di- 
agonalizc the Hamiltonian in the basis of the elements of 
the group {/, A'} to restore complex conjugation. Restor- 
ing complex conjugation amounts to letting the resulting 
5') have the property that 

k\^)=e'^\^), (13) 

where x is an arbitrary phase factor. Note that a general 
complex deformed HP determinant docs not have this 
property. 

B. Spin Projection Operators 

Due to the historical significance of spin projection, 
we here present a more detailed discussion of the form 
we use. 

We can recover eigenfunctions of by demanding that 
the projected wavefunction be invariant to the axis of 
spin quantization, as first proposed by Percus and Roten- 
berg in the context of angular momentum projectioui^ 
The form we use for the spin projection operator is simi- 
lar to the one we discussed for non-abelian groups in the 
preceding section. Explicitly, we use the operator 

P:^^^\s;m){s-k\ (14) 

= ^ / df2A^4(f^)i?(f^), (15) 

R{n) =e'°'^'e''^^ye'^^% (16) 

where -Df„j,(51) — [s] ra\R{Vt)\s\ k) is a Wigner rotation 
matrix and |s;m) is an eigenfunction of S*^ and Sz with 
eigenvalues denoted by s and m, respectively. Here, O = 
(a, /?, 7) is the set of Euler angles. 
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-^mk obeys the same rules as the transfer operators 
discussed previously: 



pS pS _ r 
mk m'k' 



ss' Sm'k Pmk' i 



{Pmkj 



(17) 
(18) 



The action of P^^j, on a deformed HF state |<i>) recov- 
ers the multi-determinantal state characterized by the 
quantum numbers s and k and writes it in terms of the 
quantum numbers s and m. As discussed above, using 
such an operator is undesirable as the wavefunction (and 
the energy) obtained depends on the value of k chosen. 
We therefore form the linear combination 



fc 



ckP:.,k\'^ 



(19) 



and determine the 2s + 1 coefficients Ck by diagonaliza- 
tion. 

We finally note that if |(f>) is an eigenfunction of Sz 
with eigenvalue m (i.e., a coUinear state), then no such 
expansion is necessary. Furthermore, the operator be- 
comes a true projection operator and can be simplified 
to 

2s -f 1 



ps 



d/3sin/3d^„(/3)e'^^«, (20) 



where dfmmiP) = (s; m|e*^'^« |s; to). In this case, the pro- 
jected wavefunctions (and their energies) depend on the 
value of {Sz) = m chosen for the deformed determinant. 
In practice, the gauge integration over the set of Euler 
angles is discretized. We have observed that the num- 
ber of grid points required to get good convergence on 
the integrals is generally small. We should note that 
the projectors we have described in the foregoing have 
frequently been used for total angular momentum pro- 
jection in nuclear physics.— It should also be observed 
that the form of the spin projector that we use is in fact 
not new in quantum chemistry. It was used by Lefebvre 
and Pralj^i*^ to provide a simpler formula for the EHF 
energy. Their work, unfortunately, remained largely un- 
noticed. 

To reduce notational clutter, we will henceforth use E 
for the energy of the projected state and will write the 
projector in a general form as 



p= dnw{n)R{n). 



(21) 



All associated eigenvalue problems are implied by our 
notation. 



C. Projected Hartree— Fock 

Given the projectors defined above, we evaluate the 
projected Hartree-Fock energy expression as 

^_ mHP\^) _ jdnw{n)mHRm<i>) ^ ^^2) 



(cf>|P|<i>) 



jdnw{n) ($|i?(i7)|$) 



Simplifying the notation by defining 

x{n) = w{n) mR{n)\<P) (23) 

and constructing the rotated Slater determinants 

mm 



\n) 



(24) 



permits us to write 

j\mx{n) {o\H\n) 



E = 



= / dny{n) {o\H\n). (25) 



Evaluating the energy thus requires evaluating overlap or 
norm matrix elements ($|i?(r2)|<I>) and Hamiltonian ma- 
trix elements {0\H\fl). Note that because R{0) ~ 1, we 
have |0) = |$). Note also that the determinants |f2) are 
defined in intermediate normalization, so that {0\fl) ~ 1. 

The norm matrix elements can be evaluated in the 
usual way for the overlap between two non-orthogonal 
Slater determinants |$) and |E:)i^ 



($|S) = detM, 



(26) 



where M is the matrix of overlaps between orbitals 
occupied in |$) and orbitals |^) occupied in |S), i.e. 



0)- 



(27) 



To evaluate the Hamiltonian matrix elements we follow 
Lowdin,— who realized that a generalized form of Wick's 
theorem holds. We have 



mm) 



Y.^^''pt? + lT.('j\\''^)pt^pf^ (28) 



where and |fc/) are the usual one-electron and anti- 
symmetrized two-electron integrals, respectively, and the 
transition density matrix elements pff" are given by 

Ptr - ^^J^ - E (kmiM-^m). (29) 

Using these formulae, we find that the overlap we need 
is given by 



($|i?(r!)|$) =detMn 
Mn = CtRnC, 

while the Hamiltonian elements are 



(30) 
(31) 



\H\n) = ^ h,kipn)k^ + \ Y.mkl){pn)k.{pn)i, 



The transition density matrices 



kl 



{<^\a\auR{n)\^) 

{mm^) 



(32) 
(33) 



5 



can be formed as 



(34) 



Here and above, C is the M x N matrix of orbital coef- 
ficients defining the occupied orbitals in |$} and Rn is 
the matrix representation of the operator i?(i7), where 
M is the number of spin-orbitals in the basis and N is 
the number of electrons in the system. We have assumed 
an ortlionormal basis for convenience and will continue 
to do so, but modifications for a non-orthonormal basis 
are straightforward. 



D. Optimization of the PHF Wave Function 

As we have written things thus far, the fundamental 
variables are the occupied orbital coefficients C. How- 
ever, because the wavefunction |$} is a single determi- 
nant, it is defined completely by its one-particle density 
matrix p, a point brought up by Lowdin in 1966,— when 
he wrote 

"Since the wavefunction Pl^") depends 
uniquely on p, the main problem in the pro- 
jected Hartree-Fock scheme is to vary p so 
that the energy becomes an absolute mini- 
mum." 

The minimization of the energy with respect to p would 
be greatly facilitated if we were able to write an explicit 
density matrix functional E[p\. Our task here is thus to 
obtain this energy functional, in a manner analogous to 
the work of Sheikh and Ring for projected Hartree-Fock- 
Bogoliuboviii Having done so, we can minimize the PHF 
energy with respect to idempotent density matrices p. 
We start by writing the density matrix as 



= CC^ 



(35) 



where we recall that C is the matrix of coefficients defin- 
ing the occupied orbitals. Let us partition this matrix in 
the form 



(36) 



where Cp is x and Cq is (Af 



A^) X N. In terms of 
these blocks of the orbital coefficients, we have 



P = 



Ppp Ppq 
Pqp Pqq 



CpCt CpCt 



c,ct 



c,ct 



(37) 
(38) 



Note that the partitioning of C is not unique, because we 
can arbitrarily mix basis functions amongst each other, 
thereby mixing rows of C. We require that Cp have non- 
vanishing determinant. We can always find a partitioning 
satisfying this requirement because we can choose to work 



in the basis of the occupied molecular orbitals of |$) for 
which Cp = 1. In practical calculations, we do exactly 
this, and the subscripts p and q denote occupied and 
virtual orbitals, respectively. 

Recall that the overlap matrix elements are given by 



($|i?(0)|$) = detMn. 
We rewrite the matrix Mo as 



Mn 



( 



(39) 



(40a) 



(Cp) 



-ip, / Cp 



(Cp)-iN,V(Ct)-\ 



Ct(Ct)-i (40b) 
(40c) 



where 



Nr 



{Ppp Ppi)^n(^p'^'^ 



(41) 



Then we can write the overlap matrix elements explicitly 
in terms of p as 



($|i?(r!)|$) 



1 



det(Nn Ppp) 



(42) 



In a similar manner, we can write the Hamiltonian 
matrix elements, computed in terms of transition density 
matrices, explicitly in terms of p. The transition density 
matrices are formed as 



t 



po = Ro (g;) M- (g^ 
Simple manipulations bring us to 



pn = Ro No (ppp Pp.) 



(43) 



(44) 



Having expressed the overlap and Hamiltonian matrix 
elements as functionals of p, we have all the ingredients 
we need to write a closed-form expression for the pro- 
jected energy in terms of p, which we recall is the density 
matrix of the deformed (unprojected) state |$). We can 
then set up a variational problem similar to that of the 
regular Hartree-Fock procedure. We simply write 



<5{i?[p]-Tr[A(p2-p)]} = 0, 



(45) 



where A is a matrix of Lagrange multipliers used to con- 
strain p to remain idempotent so that j^) remains a sin- 
gle determinant. This variational ansatz is equivalent to 
the condition that 

[JP,p]=0, (46) 
where is an effective Fock matrix given by 

d 



kl 



dpik 



E[p]. 



(47) 
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Wc provide explicit expressions for the matrix elements 
of in the appendix. 

Equation |46] can be solved by what wc may call the 
PHF eigenvalue equations 

:FC = Ce (48) 

where e is a (diagonal) matrix of orbital energies and C 
here corresponds to the full matrix of orbital coefficients 
(including virtual orbitals) . We use the eigenvectors C of 
J- to construct a new guess of the density matrix p, just 
as in standard Hartree-Fock. We remind the reader that 
these equations have been here derived in the special case 
of an orthonormal basis and must be modified slightly in 
the case of a general non-orthogonal basis. 

Thus, the optimization of the projected Hartree-Fock 
state uses the following algorithm: 

1. Construct an initial broken symmetry guess for |<f>) 
and form p. 

2. Compute the overlap matrix elements ($|i?(f2)|$) 
and form the function yi^i). 

3. Form the transition density matrices po and con- 
tract them with one- and two-electron integrals to 
evaluate the Hamiltonian matrix elements (0|i7|r2). 

4. Form the effective Fock matrix T and diagonal- 
izc it to obtain new eigenvectors and thence a new 
density matrix p. 

5. Test for convergence. If the density matrix is not 
converged, return to step 2. 

At any time the energy can be evaluated from 

= y" A^y{^) (0|if|rj). (49) 

We emphasize that the scaling of our PHF implemen- 
tation with respect to system size is the same as that 
of Hartree-Fock. At each gauge point $7, we form an 
effective one-body Hamiltonian by contracting two- 
electron integrals with transition density matrices (whose 
formation is trivial) . Having formed ^ by integrating JFq 
over the gauge angle fi, we obtain a new reference deter- 
minant by diagonalizing T . While the cost of our PHF 
procedure is moderately higher than that of Hartree- 
Fock because we must integrate over the gauge angle, this 
integration generally does not require too many points 
and is trivially parallelizable in any event. Reference 
[H establishes that the size of the gauge integration grid 
needed scales only weakly with the size of the system. 

It is important to note that the Brillouin-like condi- 
tion of Eqn. Iini defines only the occupied-virtual part of 
the effective Hamiltonian just as in standard Hartree- 
Fock theory. When working in the molecular orbital ba- 
sis as we do, this corresponds to defining T^q and T q^. 
The occupied-occupied and virtual- virtual parts of ^ are 
not defined by the Brillouin condition; accordingly, with- 
out further modifications the entire effective Hamiltonian 



vanishes at convergence. In order to cleanly separate oc- 
cupied and virtual orbitals and thereby ensure efficient 
convergence, the occupied-occupied and virtual-virtual 
parts must be defined as well. We therefore define the 
occupied-occupied and virtual-virtual parts of ^ to be 
the corresponding parts of F, the standard Fock operator 
constructed from the density matrix p. In the molecu- 
lar orbital basis, this corresponds to defining — Fpp 
and ^ qq = 'Pqq] iu a general basis, this corresponds to 
writing T ^ T + pFp (1 - p)F(l - p). This ad hoc 
choice has proven to be efficient and reliable thus far but 
is not guaranteed to work in all cases. One should note 
that the orbital energies from our PHF equations do not 
have obvious physical meaning. 

E. Nomenclature 

Before we discuss our results, let us briefly clarify the 
nomenclature we use in this manuscript. 

We will use PUHF for spin projection on the unre- 
stricted Hartree-Fock (UHF) in a projection after vari- 
ation manner. We write SUHF (equivalent to the stan- 
dard spin-projected EHF) for variation after spin projec- 
tion on a determinant which breaks 5^ symmetry and 
SGHF for variation after spin projection on a determi- 
nant which breaks both SP' and Sz symmetry. We can 
also break complex conjugation symmetry (see Ref. [ll) 
to form KUHF, KGHF, KSUHF, or KSGHF, where KS- 
GHF, for example, means that we restore complex con- 
jugation symmetry as well as and Similarly, we 
can restore point group symmetry to form, for instance, 
C2-SUHF, which restores and makes the wavefunction 
transform as one of the irreducible representations of the 
group C2. 

III. RESULTS 

Our PHF equations have been implemented in the 
Gaussian^! program package. For this pilot work, we 
will generally be interested in spin projection in small 
systems with small basis sets, though we also consider 
complex conjugation and point group symmetry break- 
ing and restoration. Larger systems with larger basis 
sets arc certainly feasible, as the method has the same 
computational scaling as Hartree-Fock. 

In order to construct an initial guess for the deformed 
determinant, we typically start with a broken-symmetry 
Hartree-Fock state. This is usually straightforward in 
the case of SUHF, but is often challenging for SGHF 
since GHF solutions different from rotated UHF are not 
common for first row molecules <^ We therefore mix the 
t-spin and 4.-spin orbitals closest to the Fermi level with 
some predefined angle small enough so as not to raise the 
energy significantly but large enough so as not to return 
to UHF. Similarly, to break complex conjugation symme- 
try, we mix the highest occupied orbital and the lowest 
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virtual orbital with complex coefficients, thereby break- 
ing complex conjugation symmetry, and we use an anal- 
ogous procedure to break point group symmetry where 
we must take care to mix orbitals belonging to different 
irreducible representations of the point group we wish 
to restore. Convergence is usually enhanced by the use 
of the Direct Inversion of the Iterative Subspace (DIIS)22, 
procedure which ideally docs not start until the energy of 
the projected state is below the energy of the symmetry- 
adapted HF state. In other words, SUHF projecting onto 
an overall singlet state should not start DIIS until the en- 
ergy is lower than the RHF energy, simply because the 
RHF wavefunction is a solution of the SUHF equations. 

Our first order of business is to verify that our SUHF 
method reproduces previous EHF results. Unfortunately, 
there are relatively few calculations available in the lit- 
erature. Nonetheless, we have verified that we repro- 
duce results from Rosenberg and Martino^^ and Klimo 
and Tihoi^. Recently, Karadakov and Cooper performed 
self-consistent projected UHF calculations based on the 
spin-coupled valence bond theoryi^ Because they worked 
within the SCVB framework, Karadakov and Cooper 
chose to keep a set of core orbitals which remained dou- 
bly occupied and symmetry adapted. The wavefunction 
is then the antisymmetrized product of the core part 
with an active part, written as the product of a set of 
non-orthogonal spatial orbitals with a single spin func- 
tion, whose coefficients are calculated variationally. In 
order to make our results directly comparable to theirs, 
wc had to carry out a constrained optimization of our 
SUHF state in a manner similar to our previous work 
on CUHF, such that symmetry breaking is only allowed 
within an active space<^i2^ We have reproduced several 
of the numbers quoted in Ref. [13 to all decimal places. 



A. Molecular Dissociation 

We start by examining the effects of projection on 
molecular dissociation curves. All calculations in this sec- 
tion use the cc-pVDZ basis set^S unless otherwise stated. 

As is well known, restricted Hartree-Fock (RHF) fails 
for the dissociation of H2, separating to an equal mix 
of two hydrogen atoms on the one hand and a hydrogen 
cation plus a hydrogen anion on the other. The ener- 
getically correct dissociation limit is recovered by break- 
ing spin and spatial symmetry in UHF. The resulting 
wavefunction is a linear combination of singlet and triplet 
and is thus spin contaminated. Restoring spin symme- 
try variationally leads to the SUHF wavefunction, which 
like UHF is energetically exact at dissociation but unlike 
UHF retains good quantum numbers along the whole po- 
tential energy curve. Near equilibrium, SUHF and UHF 
differ significantly, as can be seen in Fig. [TJ Breaking and 
restoring symmetry under K or Sz variationally improves 
the results near equilibrium even further. 

While PHF yields excellent results for the dissociation 
of H2 , it misses some of the effects of dynamical correla- 
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FIG. 1. Potential energy curve for H2 dissociation in the cc- 
pVDZ basis set. 



tion, which we illustrate by considering the dissociation 
of N2, as shown in Figs. [2]-|4l The projected HF meth- 
ods generally go to a dissociation limit slightly below 
that of UHF. Near equilibrium, breaking and restoring 
complex conjugation symmetry accounts for most of the 
correlation available in CASSCF(10/8), though KRHF 
dissociates to a limit much too high. Spin projection to 
give SUHF or SGHF is comparable to complex conjuga- 
tion projection at dissocation, but generally inferior at 
equilibrium. Combining the two in KSUHF and KSGHF 
gives energies below the (variational) CASSCF. All of 
these curves, however, are far above the coupled cluster 
singles and doubles (CCSD) based on the UHF reference, 
even with this small basis set. It is interesting to note 
that the KSGHF solution nearly parallels the UCCSD 
solution; the non-parallelity error (maximum deviation 
— minimum deviation) is 8 kcal/mol. 

A second illustration of the missing dynamical corre- 
lation is provided by considering the dissociation of H3. 
Here, we arrange the three atoms on the corners of an 
equilateral triangle and expand the triangle symmetri- 
cally. The Hartree-Fock ground state in this case is non- 
coUinear (i.e. is an eigenfunction of neither nor Sz)- 
Figure [5] shows the dissociation, projecting onto s = 1/2. 
Even in this simple three-electron system, the spin pro- 
jection clearly misses a significant portion of the total 
correlation. 

An application of spatial symmetry restoration is pro- 
vided in Fig. ini where the isotropic expansion of a square 
of hydrogen atoms is shown. Wc have performed these 
calculations with an uncontracted ST0-6G basis set, and 
wc have restored point group symmetry in the framework 
of the abclian C4 subgroup of the full Di^h symmetric 
group. In particular, the lowest energy C4-SUHF solu- 
tion corresponds to B symmetry. One can observe that 
symmetry breaking and restoration of the point group 
of the molecule provides additional correlation energy to 
SUHF, though the improvements relative to SUHF de- 
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FIG. 2. Effects of spin projection on the potential energy 
curves for N2 dissociation in the cc-pVDZ basis set, compared 
to the CASSCF(10/8) reference. 



FIG. 4. Effects of simultaneous spin and complex conjugation 
projection on the potential energy curves for N2 dissociation 
in the cc-pVDZ basis set, compared to the CASSCF(10/8) 
reference. 
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crease as the system is pulled apart. 

Unlike standard Hartree-Fock and most previous EHF 
calculations, we can calculate all the m components for 
a state with spin s. In the case of SGHF, all these states 
are degenerate, as they should be. This degeneracy is 
lost with SUHF, but a SUHF state can be found for each 
pair of quantum numbers (s, m). Figure [7] shows the dis- 
sociation of O2 to the singlet and the 771 = and 771 = ±1 
triplets, computed with a minimal (ST0-3G) basis. The 
UHF dissociation on the triplet (more precisely, on the 
777 = 1) curve is qualitatively right at equilibrium but 
has a bump before breaking spatial symmetry and dis- 
sociating from above. The lowest energy UHF solution 
at dissociation has 771 = 0. We have two different SUHF 
curves for s = 1. The 771 = 1 curve is qualitatively rea- 
sonable near equilibrium and has no bump, but goes to a 
limit above the UHF curve for the same 777. That UHF is 



below SUHF in this case is presumably because the UHF 
is contaminated by the quintet solution, which has lower 
energy at dissociation than does the triplet. The tti = 
follows the UHF "singlet" and is more correct at dissoci- 
ation. Note that while SGHF has the attractive feature 
of making all components for a given spin state degen- 
erate, finding SGHF solutions is technically demanding, 
in part because GHF solutions are generally unavailable 
but also because there are sometimes multiple SGHF so- 
lutions similar in energy. 

Another interesting aspect that has not been thor- 
oughly studied in the literature before is the basis set 
dependence of the correlation energy recovered by pro- 
jection methods. In Fig. |S1 we show the difference of 
the SUHF, KSUHF, and CASSCF(10,8) energies with 
respect to UHF as a function of the internuclear separa- 
tion with basis sets of increasing size. As it is evident 
from the figure, the correlation energy recovered is al- 
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FIG. 6. Potential energy curve for the dissociation of square 
H4 in the uncontracted STO-6G basis. C4-SUHF denotes 
the combination of spin restoration on a coUinear deformed 
determinant and restoration of C4 rotational symmetry. 
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FIG. 7. Potential energy curve for O2 dissociation in the 
ST0-3G basis set. 

most independent of the basis set size for both SUHF 
and CASSCF. This suggests that the correlations we re- 
cover are mostly static in nature. 

B. Size Consistency 

What is not so clear from the dissociation curves we 
have shown is that projected Hartree-Fock is not usu- 
ally size consistent. The magnitude of the error may be 
fairly large. For example, the KSUHF dissociation limit 
of N2 is roughly 27 mH (17 kcal/mol) higher than twice 
the KSUHF energy of the nitrogen atom, with compara- 
ble errors for KUHF and SUHF. The origin of this size 
inconsistency is simply that SUHF on N2 demands that 
the entire system be a singlet, but does not independently 
project any particular spin component onto the individ- 
ual nitrogen atoms, which should each be a quartet. 



FIG. 8. Energy relative to UHF for the dissociation of N2. 
Notice that both CAS(10/8) and SUHF show essentially no 
dependence on the size of the basis, while KSUHF has some 
small basis set sensitivity; cc-pvqz results are only shown for 
SUHF. 

In principle, a spin-projected Hartree-Fock method 
could be made size consistent with the aid of local projec- 
tion operators. That is, suppose we wish to dissociate a 
molecule to fragments A and B. Writing |$^) for 
the broken symmetry reference for fragment A (fragment 
B), we could imagine writing 

I'i') = Pab[{Pa\<^a)) ® (PbI^b))] 

, . (50) 

= PabPaPbI'^a'^b) 

where Pa and Pb are projection operators which restore 
symmetry only on the individual (noninteracting) frag- 
ments, and Pab subsequently restores symmetry on the 
whole system. We could, for example, dissociate O2 to 
two triplet oxygen atoms with the aid of Pa and Pb , and 
put the whole system into a spin triplet via Pab- This 
wavefunction is explicitly size consistent. 

While formally this looks like projection operators act- 
ing on a broken symmetry determinant to restore sym- 
metry, this approach lies outside the conventional PHF 
framework because we have relied on local projectors Pa 
and Pb in addition to a global projector Pab- Unfortu- 
nately, it is far from clear how best to define such local 
projectors except for noninteracting subsystems. The re- 
lation between the method sketched above and PHF is 
closely akin to the relation between AGP (which is not 
size consistent) and the size consistent Antisymmetrized 
Power of Strongly Orthogonal Geminals^ wherein one 
creates orthogonal subsets of orbitals and introduces 
geminals which pair orbitals only within a subset. 

C. The Curse of the Thermodynamic Limit 

We have already seen that projected HF is not size 
consistent. Neither is it size extensive. We demonstrate 
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this by examining equally spaced rings of hydrogen atoms 
with an internuclear distance of 1.8 bohr. As the number 
of atoms in the ring gets large, the ring approaches a 
periodic chain. All calculations are done in the ST0-6G 
basis set and for an overall spin singlet state. 

In Fig. [S] we show the energy per atom as a function 
of the ring size. While for small rings UHF, PUHF, and 
SUHF all differ noticeably, the three curves converge to 
the same energy per atom as the rings become large. Fig- 
ure IHl also shows that as the number of particles becomes 
large, both PUHF and SUHF give the same improvement 
in total energy over UHF. In other words, both PUHF 
and SUHF provide the same correlation energy relative 
to UHF, (about 140 mH), but the physically meaning- 
ful correlation energy per particle vanishes for large N. 
This observation was made by Mayer and Kertesa^l and 
by Castaho and Karadakov^ in a PPP model of finite 
polyene chains, and our results show the behavior they 
described. In the limit N ^ oo, where TV is the number 
of particles, restoring spin symmetry has no effect on the 
energy per particle. This is what we mean when we say 
that projected Hartree-Fock is not size extensive. 

The rationale behind the failure of projected methods 
to go beyond broken symmetry mean-field theory in the 
thermodynamic limit is as follows. For finite systems, 
the energy lowering due to projection comes about from 
the fact that the matrix elements ($|H_R|$) and ((f>|^|<f>) 
arc non-zero, where R are rotation matrices, as explained 
in Section IHl As the system size increases, those matrix 
elements tend to zero. If they are zero, then the Haniilto- 
nian becomes diagonal in the basis of the rotated states, 
and thus there is no chance for energy lowering to occur. 
In other words, the Goldstone manifold (defined in Ref. 
[l|) remains degenerate but becomes non-interacting. 



D. Singlet-Triplet Splittings 

Spin splittings (the difference in energy between the 
lowest energy states of various spins s) are important 
quantities in both organic and inorganic chemistry. They 
are connected to magnetic coupling constants in the 
Heisenberg Hamiltonian, and can be unfortunately chal- 
lenging to evaluate with simple mean-field methods. 
Density functional theory is well known to have difficul- 
ties with spin splittings, and unrestricted Hartree-Fock 
(UHF) also often fails, particularly when spin contami- 
nation is large.— 

One can attempt to remedy the failings of UHF by pro- 
jecting the UHF wavefunction onto spin eigcnfunctions 
(that is, to use a PAV approach), but in many cases this 
projected UHF is still inadequate. The question wc wish 
to address here is whether variation after projection in 
our PHF scheme offers more accurate results. We will 
focus on singlet-triplet splittings, as they are the most 
commonly studied. 

In Table U we show the computed singlet-triplet split- 
tings for a variety of small diatomic molecules as well 



as methylene (CH2)r^ trimethylmethane (TMM)^ and 
0-, m-, and p-benzyne4^ All calculations use the cc- 
pVTZ basis set, and we define the singlet-triplet split- 
tings as simply 



Ail/57^ = Ej" — Eg 



(51) 



where Et and Es arc the energies of the triplet and 
the singlet states, respectively. We have included results 
using CUHF, taken from Ref. i. In CUHF, a triplet 
state is obtained by optimizing the HF determinant sub- 
ject to the constraint that it remains a spin eigenfunc- 
tion. For singlet states, CUHF allows symmetry breaking 
only between two orbitals. PCUHF corresponds to spin- 
projection on CUHF performed with a PAV approach. 

For the small diatomic molecules, both SUHF and 
PUHF provide excellent agreement with experiment. As 
the systems become larger (but not too large; see above) 
SUHF and PUHF differ more, and SUHF is more accu- 
rate. Generally, we observe that spin projection yields 
meaningful improvements when the spin contamination 
in the UHF is large only for the singlet state, but has less 
to offer when the UHF has significant spin contamination 
in both the singlet and triplet states. This is demon- 
strated by the benzynes, where SUHF, though more 
accurate than UHF or PUHF, is still far from experi- 
ment. Complex conjugation restoration has, as one might 
expect, much smaller effects on the calculated singlet- 
triplet splittings. We should point out that dynamic cor- 
relation can have large effects on calculated singlet-triplet 
splittings which cannot be captured by SUHF and SGHF 
as these methods account primarily for static correlation. 

Note that semi-empirical corrections to calculated 
singlet-triplet splittings with spin contaminated wave- 
functions such as UHF and KUHF can be extracted as, 
for example;^"— 



AE: 



ST 



^ Ej" — Eg 



(52) 



where Et and Eg are the energies of the spin contam- 
inated triplet and singlet determinants and (S'^)t and 
{S^)s are the respective expectation values of S^. If the 
wavefunctions used are spin eigcnfunctions, Eqn. (52] re- 
duces to Eqn. 1511 While the results for spin contam- 
inated states are closer to experiment, the prescription 
is neither unique nor first principles. Applying this cor- 
rection to KUHF leads to a MAE of 10.1 kcal/mol, still 
worse than KSUHF. 



IV. CONCLUSIONS 

Symmetries and good quantum numbers are of crit- 
ical importance in many finite quantum systems. Un- 
fortunately, approximate variational solutions to the 
Schrodinger equation need not respect the same symme- 
tries as does the exact solution. Forcing them to do so 
reduces the variational ficxibility of the model, which is 
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FIG. 9. Left Panel; Energy per atom in an equally spaced hydrogen ring. Right Panel: Energetic improvement of spin 
projection over UHF. Both panels use one over the number of atoms for the abscissa; the thermodynamic limit is — ^ 0. 
Color online. 



TABLE I. Singlet triplet splittings (kcal/mol) from several methods. "TMM" refers to trimethylmethane, while ME is the 
mean error (theory — experiment) and MAE is the mean absolute error. 



Molecule 


UHF 


KUHF 


CUHF 


PUHF 


PCUHF 


SUHF 


KSUHF 


Expt 


NH 


19.4 


18.6 


21.0 


38.2 


41.5 


33.6 


31.6 


39.0 


OH+ 


25.9 


25.0 


27.4 


50.9 


54.1 


45.8 


43.4 


50.6 


O2 


15.8 


14.6 


16.1 


26.8 


30.6 


20.6 


24.2 


22.6 


NF 


19.7 


18.6 


20.6 


38.1 


40.8 


32.3 


31.0 


34.3 


CH2 


16.9 


15.9 


15.4 


18.9 


15.7 


15.6 


14.0 


9.4 


TMM 


23.7 


15.6 


7.5 


28.6 


10.5 


19.1 


18.0 


17.7 


o-benzyne 


-15.8 


-12.5 


-9.8 


-42.6 


-24.5 


-51.4 


-48.7 


-38.0 


m-benzyne 


28.3 


13.0 


12.8 


45.0 


12.8 


2.2 


-9.1 


-20.6 


p-benzyne 


-10.1 


-5.2 


0.1 


-24.0 


-0.5 


-28.2 


-22.9 


-3.5 


ME 


1.3 


-1.0 


0.0 


7.6 


7.7 


-2.5 


-3.4 




MAE 


17.5 


15.5 


15.9 


13.4 


9.3 


9.2 


7.4 





not ideal. On the other hand, symmetry breaking from 
mean-field theories generally indicates the failure of the 
mean-field approximation and the emergence of strong 
correlation. One can take advantage of this symmetry 
breaking to obtain energetically reasonable broken sym- 
metry mean-field wavefunctions, though the wavcfunc- 
tions obtained are of poor quality. By using projection 
operators, one can restore the broken symmetries and 
recover the associated quantum numbers while obtaining 
strongly correlated wavefunctions. Moreover, this can be 
accomplished without leaving the independent particle 
picture that allows one to easily grasp the physics in the 
wavefunction. 

While sclf-consistently restoring the symmetry of 
mean-field wavefunctions is formally attractive, the prac- 
tical realization of it is not trivial. On the other hand, 
a projection-after-variation approach is simple, but may 
lead to severe problems with the projected wavefunction. 
A self-consistent approach also allows one to break sym- 



metry in a deliberate way, that is, one can find a broken 
symmetry reference determinant even in cases when the 
standard HF method leads to a symmetry-adapted state. 

Our earlier work showed how to optimize the energy of 
the projected wavefunction when the broken symmetry 
reference state was a Hartree-Fock-Bogoliubov determi- 
nant that broke particle number symmetry. Here, we 
do the same when the reference determinant conserves 
electron number. There are two reasons to present this 
new set of PHF equations. First, due to the historical 
significance of the PHF problem in quantum chemistry. 
Second, the PQT equations as presented before lead to 
indeterminancies for systems whose deformed reference 
state has exact zero occupations, as is the case when 
particle number symmetry is not broken. The pairing in- 
teraction in PQT assigns a non-zero occupation to every 
natural orbital, but in practice the occupation numbers 
may become sufficiently small in a large enough basis 
that the equations cannot be accurately solved in double 
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precision arithmetic. The ideas presented here can be 
extended to remedy this deficiency. 

The PHF formulation here presented is easy to 
code and numerically robust. The resulting projected 
Hartree-Fock scheme is conceptually simple and compu- 
tationally affordable, and generally leads to good qual- 
ity multirefcrence wavcf unctions. Size consistency and 
size extcnsivity remain problems, but otherwise projected 
Hartree-Fock offers a relatively black box treatment of 
strong correlation while staying within the mean-field 
picture without sacrificing good quantum numbers. 



The derivation of the matrices Yq and JF is straight- 
forward but cumbersome. Our final expression for Yn 
is 

Yn^yn- J dn'y{Q')yn' (A7) 



where 



yn - I I No + No [p.,p ppj Tin- (A8) 
\Pgp / 

If we define Fq = h + G^, we can write as 
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Appendix A: The PHF Effective Hamiltonian 

Here, we provide the explicit expressions required to 
evaluate the matrix elements of the effective Fock matrix. 
The PHF energy expression can be conveniently written 
as 



E = J dny{n) I^TT[h pn] + ^TiiGn pn] 



(Al) 



where h is the matrix of one-electron integrals and Go 
is 



{GnU = Y.('j\\kl)ipn) 



(A2) 



The effective Fock matrix is defined in terms of deriva- 
tives of the PHF energy expression with respect to matrix 
elements of the density matrix (of the deformed state) as 



J'ki — 



d 



dp. 



Ik 



-E[p]. 



^ is Hermitian by construction, that is. 



(A3) 



(A4) 



We follow Sheikh and Ringi^ in writing the derivatives 
of the function yi^C) in the form 



d 



dpki 



yin) = y{Q) (Yo)ife. 



(A5) 



Using this, we can write the effective Fock matrix as 



dny(n){(Ynhj 



Tr(/ipo) + ^Tr(GoPo) 



n)ik 



kl 



{pQ.)kl^- 



(A6) 



(A9) 



where 



J^n = ^YoTr[(/i + Fo)po] 



Nj2 [ppp Ppq) Ff2(l - po)R.n 



(1 - po)FoRo 



Nr 



(AlO) 



Note that since No is an x iV matrix, both y^ and 
parts of J^Q, are given as a sum of an M x A'' and an x 
M matrix. These matrices are correctly understood as 
M X M, padded to the right or below by zeroes. The fact 
that {Tnjqq and {ynjqq vanish is simply a consequence 
of the fact that the energy is independent of pqq. 

In the molecular orbital basis of |$) in which we work 
in practice, we have 



P = 



1 0' 

,0 0, 



(All) 



The subscripts p and q here will be replaced by o and 
V to emphasize that they denote occupied and virtual 
orbitals, respectively. Writing 



o — 



1 1^7; o R" 



(A12) 



where we have suppressed the dependence on for 
brevity, we obtain 



No = 



(A13) 



where Rq^^ is the inverse of Roo- The transition density 
matrices are then 



Pn 



1 0' 

RuoRoo , 



One finds that 

yn 



^00 "^ov 



(A14) 



(A15) 
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Finally, we have 

J=a = ^YoTr \{h + Fo)po] + (A16) 

where the individual blocks of are 

(^o)oo = 0, (A17a) 
[Rt,t, — Rt,oRoo'"Roii] , (A17b) 

(■^o)uo = F„o — R^oRqq"'"Foo (A17c) 

= 0. (A17d) 



Using Eqn. IA7l and the fact that 2/(0) integrates to unity, 
it is clear that the occupied-occupied and virtual-virtual 
blocks of Yo vanish, and thus so do these blocks of T^. 
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